Zero-one survival behavior of cyclically competing species 



Maximilian Berr 1 ' 2 , Tobias Reichenbach 1 , Martin Schottenloher 2 , and Erwin Frey 1 
1 Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), 
LMU Munchen, Theresienstrafie 37, 80333 Munchen, Germany 
2 Mathematisches Institut der LMU Munchen, Theresienstrafie 39, 80333 Munchen, Germany 

(Dated: January 30, 2009) 

Coexistence of competing species is, due to unavoidable fluctuations, always transient. In this 
Letter, we investigate the ultimate survival probabilities characterizing different species in cyclic 
competition. We show that they often obey a surprisingly simple, though non-trivial behavior. 
Within a model where coexistence is neutrally stable, we demonstrate a robust zero-one law: When 
the interactions between the three species are (generically) asymmetric, the 'weakest' species survives 
at a probability that tends to one for large population sizes, while the other two are guaranteed to 
extinct. We rationalize our findings from stochastic simulations by an analytic approach. 
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Ecological systems are composed of a large number 
of different interacting species pQ. Competition between 
them basically affects the probability of individuals' re- 
production as well as death. However, such birth and 
death processes also possess a considerable degree of 
stochasticity, which induces fluctuations that ultimately 
result in species extinction [2j[3]. In this respect, further 
understanding of the conditions and mechanisms that en- 
able the huge observed biodiversity is subject of a large 
body of work in contemporary theoretical ecology and 
biological physics. It involves the challenging problems 
of characterizing out-of-equilibrium steady states in the 
presence of intrinsic fluctuations [3] , nonlinearities [5] as 
well as nontrivial interaction networks [6] . Recent exper- 
iments on colicinogenic microbes have proven the impor- 
tance of cyclic, 'rock-paper-scissors' like, competition [7]. 
Such cyclic dynamics also governs, e.g., certain lizard 
populations [8] and coral reef invertebrates [9J. Theo- 
retical studies have mainly focussed on identifying con- 
ditions under which cyclic competition leads to main- 
tained diversity, employing, e.g., a time-scale framework 
to distinguish stable from unstable coexistence [10]. A 
supporting role of self-forming spatial patterns has been 
underlined generally [TQl [TTJ [12] , although in certain sit- 
uations it may not be necessary [13] or occasionally even 
harmful [13] . 

In contrast, little is known about the fingerprints of 
extinction. E.g., in the E.coli experiments [7], when the 
population is well-mixed, a strain which is resistant to the 
poison without producing toxin itself remains as the only 
survivor after a short transient. Is this behavior robust? 
If so, why does one species reliably outcompete the other 
two although all three species together display a cyclic 
dynamics, where each out competes another but is itself 
beaten by the remaining one? What is the influence of 
unavoidable fluctuations? 

In this Letter, we approach these ecologically impor- 
tant and physically insightful questions by investigating 
cyclic competition of three interacting species, referred to 



as A, B, and C. Aiming at a broad and general applica- 
bility of our model and results, we consider the following 
simplified, paradigmatic interactions: 

A + B A + A, 
B + C ^ B + B, 
C + A^C + C. (1) 

Hereby, A outperforms B at a rate fc^, while B beats 
C which outcompetes A in turn, at rates ks resp. kc- 
Recently, it has been shown that a population of N such 
interacting individuals eventually ends up in one of the 
(absorbing) states where only one species survives [15] 
US] . The mean time T for extinction is proportional to 
the system size N, T ~ N, indicating that extinction is 
solely driven by fluctuations p~5j [16]. Therefore, which 
one of the species survives is subject to a random process. 
If the interaction rates are equal, i.e., the three species 
are symmetric, all have an equal chance of surviving. 

Here, we investigate the generic case where the com- 
peting species do not obey a symmetry. To gain intuition 
for the system's behavior, we discuss predictions by the 
rate equations (RE) first. The RE describe the deter- 
ministic time-evolution of the densities a, b and c of the 
three species, as may arise when fluctuations are negli- 
gible, e.g., when the population size N is large. For the 
reactions ([!]), the RE are given by 

a = a(kAb — kcc) , 
b = b(kBC — kAd) , 

c = c(kca — ksb) . (2) 

They conserve the number N of interacting individuals: 
the densities fullfill the relation = 1, spanning the 

simplex 63 as the phase space (see Fig. [I]). Its corners 
represent three absorbing fixed points, where only one 
species remains. In addition, Eqs. ([2| possess a reactive 
fixed point F, located at (a*, b* , c*) = (fc#, fees &a)/(&U + 
ks + kc) which corresponds to species coexistence. In the 
following, we use the time-scale normalization condition 
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FIG. 1: The phase space S3. We show the reactive fixed point 
F, the center Z, as well as a stochastic trajectory (red/light 
gray). It eventually deviates from the 'outermost' determin- 
istic orbit (black) and reaches the absorbing boundary. For 
the distances Aa, As and Ac (blue/dark gray) see text. Pa- 
rameters are (/ca, kc) = (0.2, 0.4, 0.4) and N — 36. 



ks + kc + kA = 1. Therewith, the parameter space, 
spanned by the rates kA, kc, also adopts the form of 
the simplex S3, see the lower parts of Fig. [2j 

The RE ([2| predict neutral stability of species coexis- 
tence, as they obey the following constant of motion: 

R = a kB b kc c kA , (3) 

which does not change in the course of the determin- 
istic time-evolution. Similar to the energy in classical 
mechanics, R singles out closed orbits surrounding the 
coexistence fixed point F, see Fig. [I] These orbits, as 
well as F, are neutrally stable to fluctuations; stochastic 
trajectories follow the cyclic behavior of the determin- 
istic orbits to a certain degree, while at the same time 
performing a random walk between them. Eventually 
they reach the boundary of the phase space and are then 
driven to one of the absorbing fixed points (c.f. Fig. [I]). 

We have performed extensive computer simulations to 
determine the influence of the reaction rates kA,ks,kc 
as well as the system size N on the probabilities P SU rv 
for each species to survive. To this end we have evolved 
the system, with initial coexistence, until extinction of 
two species occurred; the average outcome over many 
such runs defined the survival probabilities. Typically, 
the system has initially been in a state corresponding to 
the center of the simplex. However, altering this starting 
point is not relevant for the results, as we show below. 

What is the influence of the population size N on the 
survival behavior? To answer this question, firstly, we 
consider the smallest population where all three species 
can 'coexist', namely N = 3. In this case, P SU rv de- 
pend linearly on the reaction rates, see Fig. [2] (a). For 
such a small system, the master equation describing the 
stochastic processes ([I]) can be solved exactly. Only the 



(a) N = 3 (b) N= 10,000 








FIG. 2: Survival probabilities for A (blue/dark gray), B 
(green/medium gray) and C (red/light gray), obtained from 
stochastic simulations as averages over 10,000 samples, (a), 
A small system, N = 3, leads to a 'law of stay-out' and linear 
dependences of P SU rv on the reaction rates, (b), Large pop- 
ulations, here for N = 10, 000, are governed by a 'law of the 
weakest' and obey, in the limit N — > 00, a zero-one behavior. 



state where one individuals of each species is present cor- 
responds to coexistence, the other states lie on the ab- 
sorbing boundary. A single process therefore immedi- 
ately leads to extinction and determines which species 
survives. The resulting survival probability for species A 
reads P^ rv = &b, the others follow analogously. A 'law 
of stay-out' arises: The species that is least frequently 
engaged in interactions (for species A, interactions occur 
at rates kA and kc) has the highest chance to survive. 
In contrast to what emerges for large populations, see 
below, this law is not strict. If ks denotes the largest 
of the interaction rates, species A is not guaranteed to 
survive, but possesses the highest probability. 

Large populations, about N > 100, display a con- 
trasting 'law of the weakest' which determines the sur- 
viving species [17]. Namely, for reaction rates fullfilling 
kA < kB,kc, species A has the highest probability of 
surviving, although A may be considered as the 'weakest' 
species: Its reproduction occurs at rate fc^ and is thus 
the slowest of the three competing species. This non- 
trivial law has previously been described, as a non-strict 
one, in Ref. [17]. Here we show that, surprisingly, this 
law becomes strict in the limit of large population sizes, 
TV -> 00 (see Fig. [2] (b) for the situation N = 10,000). 
In this limit, P^ rv — ► 1 in the region kA < ks,kc, and 
P^rv ~~ > otherwise; a 'zero-one' behavior arises. Three 
regions emerge in parameter space, depicted in the lower 
part of Fig. [2] (b). In each of them one distinct species is 
guaranteed to survive, while the others go extinct. 

The transition from the 'law of stay-out' to the 'law 
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FIG. 3: The limit of large populations. We show the survival 
probabilities of A (open symbols) and B (filled) depending 
on ks, for kc = 0.35 and different system sizes: N = 24 (red 
squares), 100 (blue circles), and 10,000 (black triangles). A 
sharpened transition emerges at k* B = 0.325 where the A- and 
B-dominated regions meet. Our analytical prediction, Eq. |7| 
is shown for p = 0.7 and N = 10, 000 (dashed line). 



of the weakest' happens gradually upon increasing the 
system size. From our stochastic simulations, small pop- 
ulations, around N < 20, are governed by the 'law of 
stay-out' (although strictly valid only for N = 3). Inter- 
mediate scenarios emerge for medium populations, about 
20 < N < 100, while survival in large ones, N > 100, is 
predominantly determined by the 'law of the weakest'. 

We have confirmed that the 'law of the weakest' be- 
comes strict for N — > oo by computing the survival prob- 
ability for different population sizes N. For illustration, 
let us focus on a one-dimensional section through the 
two-dimensional parameter space, connecting a region 
where A is 'weakest' to another one where B is 'weak- 
est'. Results are shown in Fig. [3j they demonstrate that 
the transition between A- and ^-survival, at the 'criti- 
cal' parameter value where both regions meet, becomes 
sharp when increasing the population size, resulting in a 
discontinuous transition. 

Further evidence for a sharp transition in the survival 
probabilities, i.e., a zero-one behavior, comes from an 
analytical approach which we have developed. It de- 
scribes the survival probabilities for large systems, where 
only the last part of the stochastic time-evolution, shortly 
before extinction, is relevant. Indeed, in large systems, 
the stochastic trajectories exhibit only small fluctuations, 
and closely follow the deterministic orbits, performing 
many turns around the reactive fixed point F (see Fig. [I]). 
Eventually, they deviate from an 'outermost' determinis- 
tic orbit, shown in black in Fig. [I] and hit the absorbing 
boundary. Then, the system is driven to one of the ab- 
sorbing states. Which one is reached depends on which 
edge of phase space the trajectory had reached before; 
each edge leads to one distinct uniform state. 

For asymmetric reaction rates, the fixed point is shifted 
from the center Z of the phase space towards one of the 



three edges, i.e., into one of the three domains shown in 
Fig. [I] All deterministic surrounding orbits are changed 
in the same way, squeezing in the direction of one edge. 
Intuitively, the absorbing state which is reached from this 
edge has the highest probability of being hit, as the dis- 
tance from the 'outermost' deterministic orbit towards 
this edge is shortest. Indeed, this behavior has been val- 
idated by the above presented stochastic simulations. 

Let us formalize the above considerations. We define 
the outermost deterministic orbit as the one orbit that 
is only a distance of 1/iV, i.e., one discrete, elementary 
step, apart from its closest edge. The distance of this out- 
ermost orbit to the edge that induces survival of species 
A is termed A^; the distances A# and Ac are defined 
analogously. Now, in the parameter region fc^ < 
where A has the highest survival probability, the distance 
A a is smallest, and therefore A a = 1/N. The other two 
distances can be obtained via the conserved quantity R, 
Eq. ([3]). For this purpose, in the following, we consider 
the (most interesting) situation where the differences be- 
tween the reaction rates &u,&b,&c are small. The out- 
ermost orbit then runs through the point c = Xa = 1/N 
and a « b « 1/2, yielding its constant of motion 



R o.O = 



1 



J\fk A 2 k B+k c 



(4) 



If we perform the same calculation at the point where the 
outermost orbit lies closest to the edge leading to survival 
of species B resp. C, we obtain 



R o.o = x k B 1 and R o.o = x k c 



resp.. Equating these expressions yields 



(5) 



X B = 2~^~ xN~^ and X c = 2~^~ xN~^ . (6) 

Most important for our purpose is the scaling of the 
distances A#,Ac in the population size N. Residing 
within the regime fc^ < kB,kc, we notice from Eqs. (|6| 
that both Xa and Xb decrease slower than 1/N. Conse- 
quently, the number of discrete steps that separate the 
outermost orbit from these two edges, given by NXb resp. 
TVAc, tend, for large populations, to infinity. Note that 
the same does not apply to A^, which we keep fixed at 
1/N. Below, we show how this scaling leads to the zero- 
one behavior of the survival probabilities. 

The probability for deviating from the outermost orbit 
and performing an elementary step towards the absorb- 
ing boundary of the phase space is, for small differences 
in the reaction rates, approximately constant along the 
orbit, let us denote it by p < 1. Now, the probability of 
leaving the outermost orbit and reaching the edge lead- 
ing to survival of species B is given by the probability 
of NXb such subsequent elementary steps, and therefore 
reads p NXs . Analogously, we obtain p NXA and p NXc as 
the probabilities for reaching the edges connected to A 
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and C, resp.. Consequently, the survival probabilities of 
the different species are given by the corresponding nor- 
malized probabilities: 
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(7) 



6 follow analogously. The above found scal- 
N\b, N\c — > oo upon N — > oo, while NXa = 1, 



surv ajl± ^ L J surv 
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imply that for large populations 
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0? P S urv ~~ * 0? an d therefore the zero-one behavior 
emerges. Note that these results have been derived from 
assuming kA < ks^kc'-, different species survive in the 
other regions of the parameter space. The overall sig- 
moidal form of Eq. ([7| and the associated width agrees 
well with numerical findings for large systems, see Fig. [3] 
However, deviations do occur in the more detailed shape 
of the survival probability. We attribute them to the 
approximations we made when deriving Eq. (|7|), namely 
that we have treated the escape step probability p along 
the outermost orbit as well as the latter's distances to 
the different edges as constant. 

Symmetries between species alter the survival proba- 
bilities. If all interaction rates are identical, all species 
clearly have the same chance of 1/3 to survive; and if two 
rates coincide, the corresponding two species can both 
have chance 1/2 of remaining. 

The above analysis suggests that the survival probabil- 
ities do not depend on the starting point, as long as the 
latter is not too close to the absorbing boundary. Indeed, 
extinction occurs due to deviations from the outermost 
orbit, any initial state therein will induce the same behav- 
ior. This expectation has been confirmed by simulations. 

The dependence of the extinction behavior solely on 
temporally late deviations from the outermost orbit is 
reminiscent of 'tail-events' in probability theory [18] 
which induce the celebrated zero-one-law originally due 
to Kolmogorov [18]. While this law cannot be directly 
applied to the present situation, mainly due to the finite 
number of steps in each of the trajectories discussed here, 
further investigations along these lines seem promising to 
deepen our understanding of zero-one behaviors. 

We have derived the zero-one behavior of survival 
probabilities, accompanied by a strict 'law of the weak- 
est', for a cyclic population model that, deterministically, 
exhibits neutrally stable coexistence. However, the above 
analysis based on scaling arguments allows us to imme- 
diately generalize the obtained results to the case where 
coexistence is (deterministically) stable. Then, stochas- 
tic trajectories are attracted to the reactive fixed point, 
with rare large deviations. Again, extinction is deter- 
mined by the behavior of trajectories close to the absorb- 
ing boundary, such that an analogous analysis as derived 
above holds. However, the latter is much harder to test 
numerically. Deterministically stable coexistence induces 
a mean extinction time which increases exponentially in 



the system size [T9] , such that computation of survival 
probabilities in large systems is hardly feasible. 

We conclude by relating our results to the E. coli exper- 
iments [7] mentioned in the beginning. Identifying A with 
the sensitive, B with the resistant, and C with the coli- 
cinogenic strain, we uncover the relation kc ^> > fc# 
[C can kill A (fast) and reproduce, while A and B can 
only reproduce if a neighboring bacterium dies (slow). 
kA and ks are then proportional to the reproduction rate 
differences of A and B resp. of B and C. The measured 
data [7] leads to kA > &b ]• The resistant strain B is thus 
'weakest' and, according to the 'law of the weakest', sur- 
vives, in agreement with experimental observations [7]. 
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